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ABSTRACT 

Tidal stellar disruptions have traditionally been discussed as a probe of the single, massive black 
holes (MBHs) that are dormant in the nuclei of galaxies. In Chen et al. (2009), we used numerical 
scattering experiments to show that three-body interactions between bound stars in a stellar cusp 
and a non-evolving "hard" MBH binary will also produce a burst of tidal disruptions, caused by a 
combination of the secular "Kozai effect" and by close resonant encounters with the secondary hole. 
Here we derive basic analytical scalings of the stellar disruption rates with the system parameters, 
assess the relative importance of the Kozai and resonant encounter mechanisms as a function of time, 
discuss the impact of general relativistic (GR) and extended stellar cusp effects, and develop a hybrid 
model to sclf-consistently follow the shrinking of an MBH binary in a stellar background, including 
slingshot ejections and tidal disruptions. In the case of a fiducial binary with primary hole mass 
Mi = 10 7 Mq and mass ratio q = M2/M1 = 1/81, embedded in an isothermal cusp, we derive a 
stellar disruption rate A* ~ 0.2 yr" 1 lasting ~3x 10 5 yr. This rate is 3 orders of magnitude larger 
than the corresponding value for a single MBH fed by two-body relaxation, confirming our previous 
findings. For q <C 0.01, the Kozai/chaotic effect could be quenched due to GR/cusp effects by an 
order of magnitude, but even in this case the stellar-disruption rate is still two orders of magnitude 
larger than that given by standard relaxation processes around a single MBH. Our results suggest 
that > 10% of the tidal-disruption events may originate in MBH binaries. 
Subject headings: black hole physics - methods: numerical - stellar dynamics 



1. INTRODUCTION 

Stars that wander too close to the MBHs that reside 
at the center of galaxies arc shredded by the tidal gravi- 
tational field of the hole. After a tidal disruption event, 
about half of the debris are spewed into eccentric bound 
orbits and fall back onto the hole, giving rise to a bright 
UV/X-ray outburst that may last for a few years (e.g. 
iReesI 119881) . "Tidal flares" from MBHs may have been 
observed in several near by inactive galaxies (|Komossa I 
l200aiEsquei et al.ll2007l ). The inferred stellar disruption 
frequency is ~ 10~ 5 yr " 1 per galaxy (wit h an order of 
magnitude uncertainty. iDonlev et al.l 120021 ) . comparable 
to the theoretical ex pectations for single M BHs fed by 
two-body relaxation (jWang fe Merritd [20041 ). 

Yet, MBHs are not expected to grow in isolation. Ac- 
cording to the standard paradigm of structure forma- 
tion in the universe, galaxies merge frequently during 
the assembly of their dark matter halos. As MBHs 
become incorporated into larger and larger halos, they 
sink to the center of the more massive progenitor ow- 
ing to dynamical friction from distant stars, and form 
bound binaries (MBHBs). In a purely stellar back- 
ground, as the binary separation decays, the effective- 
ness of dynamical friction slowly declines, and the pair 
then "hardens" via three-body interactions, i.e. by 
capturing stars that pass close to the holes and eject- 
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ing th em at much hi gher velocities (e.g. iBegelman et al 
19801 : lQuinlaiJH996t IVoTonteri et al. Il2003t ISesana~ 



et al 



20061 ). If the hardening continues sufficiently far, pos- 



sibly driven by efficien t stellar relaxation pro cesses in a 
triaxial potential (e.g. IMerritt fe Poon| 12004) or i n the 
presence of massive pert urbers (e.g. IPerets et al.l [20071 : 
IPerets fe Alexande r 2008) , or by dissipative gaseous pro 
cesses (e.g. IColpi fc Dottill2l)09f l. gravitational radiation 



losses finally take over, and the two MBHs will coalesce 
in less than a Hubble time (e.g . IMerritt fc Milosavlievic I 
l200llSesana et alJ2505Ll200l . In lChen et al.l (120091) . we 
used scattering experiments to show that gravitational 
slingshot interactions between a non-evolving, unequal- 
mass hard binary and a bound stellar cusp will inevitably 
be accompanied by a burst of stel lar tidal disr u ptions . 
Our work differed from those by llvanov et al.l (|2005l) . 
who developed an analytical theory of the secular evolu- 
tion of stellar orbits in the gravitational field of a MBHB, 
and by IChen et al.l ([20081 ) , who argued that stellar dis- 
ruption rates by MBHBs fed by two-body relaxation 
would be smaller than those expected for single MBHs. 
Our numerical experiments revealed that a significant 
fraction of stars initially bound to the primary hole are 
scattered into its tidal disruption loss cone by resonant 
interactions with the secondary hole, close encounters 
that change the stellar orbital parameters in a chaotic 
way. 

In this paper we continue our investigations of stellar 
disruptions by MBHBs embedded in bound stellar cusps. 
We develop a hybrid model that self-consistently follows 
over time the shrinking of an MBH binary, the evolu- 
tion of the stellar cusp, and the stellar disruption rate. 
The plan is as follows. In § [21 we introduce the basic 
theory of stellar disruption processes by MBHB systems. 
We describe our numerical scattering experiments in § [3l 
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and discuss our results for different binary parameters as 
well as the effect of general relativistic corrections in § |4j 
A detailed study of the properties of disrupted stars is 
carried out in § [5] As a first step towards understanding 
the dependence of stellar consumptions on the parame- 
ters of the system, in §|6]wc fix the binary scmimajor axis 
and its eccentricity, and calculate the stellar disruption 
rate in the stationary case. In §0 we present our hybrid 
model and calculate the disruption rates for an evolving, 
shrinking MBHB. Finally, we summarize and discuss our 
results in § HI 

2. BASIC THEORY OF STELLAR DISRUPTIONS 

Consider an isotropic background of stars all of mass 
to* and radius r», centered on an MBH. Let fy(r) be 
the total gravitational potential at radius r, and ft = 
^(Mbh/to*) 1 ' 3 the tidal disruption radius, 
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(1) 



The phase-space region of specific energy E* and specific 
angular momentum J* bounded by 



J 2 c (E*,r t ) = 2r 2 [E*~y(r t )} 



(2) 



is populated by stars on orbits crossing rt, and thus sus- 
ceptible to tidal disruption. We name this cone-like re- 
gion of phase space the "tidal loss cone" . Whether the 
tidal loss cone can be emptied by stellar disruption de- 
pends on the efficiency of stellar relaxation. Let T r (E*) 
be the relaxation timescale of stars with specific energy 
E*, P*(_E*) their orbital period, and J C {E*) the spe- 
cific angular momentum of a circular orbit with energ y 
E*. In the "pinhole limit" (jLightman fc Shapiro! 119771 ). 
P*(E*)/T r (E*) > J? c (E*,r t )/J 2 (E*), a star can ran- 
dom walk in and out of the tidal loss cone within one 
orbital period, and the tidal loss cone remains almost 
full despite tidal disruptions. In the "diffusion limit" , 

P*(£*) /T r (E*) < J&(E*,r t )/J2( E *)> the tidal loss cone 
is emptied after a single orbital period, and stars diffuse 
into the loss cone on the relaxation timescale. Assume 
now that the central primary hole of mass M\ forms a 
binary pair with a secondary hole of mass M2 < M\ , and 
let a be the semimajor axis of the system. The E* — J* 
region of phase space bounded by 



J? c (E*,a) = 2a 2 [E* - *(o)] 



(3) 



is composed of orbits that are either inside or intersect 
a sphere of radius a. If the binary is "hard", a star on 
such orbit will undergo a three-body interaction with the 
MBHB, so we refer to the phase space defined by equa- 
tion Q as the "interaction loss cone". Three-body in- 
teractions perturb the energy and angular momentum of 
"intruder" stars, acting as an additional source of stellar 
relaxation. If three-body relaxation occurs in the dif- 
fusion regime, the stellar consumption rate will be en- 
hanced. 

To proceed further, we must first define some char- 
acteristic scales of a MBHB system. Recent numer- 
ical simulations have shown that three-body interac- 
tions between the binary and intruder stars result in 
significant energy exchange when the total stellar mass 
within the binary orbit is comparable to or smaller than 



the mass of the secondary hole (|Baumgardt et all 120061 : 
iMatsubavashi et all 120071) . We denote with ao such a 
critical binary separation: the binary shrinks by dynam- 
ical friction when a > ao, and by three-body p r ocesse s 
at smaller separations. Following ISesana et al.l (|2008[) . 
we assume that the stellar distribution follows a double 
power- law with break radius ro , defined as the radius of 
the "sphere of influence" containing a mass in stars equal 
to 2M\. For r > ro, the stellar density profile follows an 
isothermal distribution, 



P*(r) 



2nGr 2 ' 



(4) 



where er* is the 1-D velocity dispersion, while for r < tq 
p* (r) on r 1 . It is easy to derive then 

ro = (3 - ~/)GM 1 /a 2 ~ 4.6 pc (3 - j)M 7 a^ 2 (5) 

and a = q^^'^ro, where M 7 = Mi/10 7 Mq, 
q = M2/M1 is the binary mass ratio, and <7ioo = 
<7*/100 km s^ 1 . Notice that the "three-body radius" 
ao is larger than t he conventiona l "hardening" radius 
a h = GAf 2 /(4cr 2 ) (|Quinlanlll996l) . The ratio between 
the tidal radius of the primary hole, r t i, and ao, 
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(6) 



where p = 1/(3 — 7), indicates that the interaction loss 
cone of a binary is much larger than the tidal loss cone 
of a single MBH. Therefore, the transfer of only a small 
fraction of interacting stars into the tidal loss cone will 
cause a large enhancement of the stellar disruption rate. 

When does the presence of a binary begin affecting the 
stellar disruption rate? Let us assume that, before the 
intrusion of the secondary hole, stellar relaxation is dom- 
inated by two-body interactions. Stars in the diffusion 
limit are bound to the primary hole, and their specific 
energy E* is related to the orbital semimajor axis a* by 
E* = —GMi/(2a*). The boundary between the pinhole 
and diffusion limits is then dictated by the condition 



P*(a*)/T r (a*) = Ji c (a*,r t i)/Jc(a,*) = ni/a* 



(7) 



Substituting into the above equation the two-body relax- 
ation timescale, 



T r (r): 
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(where In A is the Coulomb logarithm) and the Keplc- 
rian orbital period of the star orbiting Mi, P*(a*) = 
2n[al/(GMi)] 1 ^ 2 , and assuming r* = R Q and to* = M , 
we can write the critical radius a c marking the boundary 
between pinhole and diffusion regimes as 



ro 



0.33 ' (3-7) 1 M 7 ' cr 100 



In A 
Tff 



-2/s 



(9) 



where s = 5 — 2j. If a secondary hole is now added to 
the system, and the interaction loss cone is not empty, 
a significant enhancement of stellar disruptions occurs 
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when the binary separation shrinks to a ~ a c . For an 
isothermal density profile and a primary hole satisfyin g 
the Mbh — c* relation, M7 = irf 00 dTremaine et al.ll2003) . 

equation © implies a c > a as long as q < 0.1 A/^ 3 , i.e. 
for unequal-mass binaries the enhancement of stellar dis- 
ruptions starts during the dynamical friction early phases 
of the binary orbital evolution. N-body simulations have 
shown, however, that the secondary hole decays from a c 
to ao and enters the three-body interaction regime on a 
timcscalc < 10 5 yr. Here, we ignore the early dynamical 
friction phases and focus on stellar disruptions induced 
by three-body scattering events. At binary separation 
a = ao, the interaction loss cone contains stars that can 
be bound or unbound to the primary. A bound star 
can interact with the binary multiple times before leav- 
ing the system, significantly increasing its probability of 
being tidally disrupted. For equal-mass binaries, the ra- 
dius of influence ro is comparable to the three-body ra- 
dius ao, and most scattering events involve stars that are 
unbound (or marginally bound). The impact of bound 
stars is more important for unequal-mass MBHBs, and 
these systems will be the main focus of this paper. 

A bound star with semimajor axis a* < a/2 never 
crosses the orbit of the secondary hole and undergoes 
a secular evolution in which its orbital eccentricity is 
excited and oscillates periodically, t he so-called "Kozai 
effect" (IKozai I U96l (Lidov I 11962 llvanov et al.l [20051: 
iGualandris fc Merritl I2009T ). The period of oscillation 
( "Kozai timescale" ) is 

2 /a \ -3/2 
TM = — (-) P(a), (10) 

(jlnnanen et alJll997t Ikiseleva et alJl!998f >. where 
P(a) = 2™ 3 / 2 [G(M X + M 2 )P 1/2 

^10 3 yr(l + a)-^A/-V^^^ /2 

is the orbital period of the binary. Since Tk (a) *C T r {a), 
the Kozai mechanism is much more efficient than two- 
body interactions at repopulating the tidal loss cone. 
However, when g < 1, r > % and the majority of 
bound stars have close encounters with the secondary 
hole that change the orbital elements of the star in a 
complicated chaotic way. In this regime: 1) numerical 
simulations are needed to give reasonable estimates of 
the tidal disruption rates; and 2) the contribution to the 
gravitational potential by background stars as well as 
stellar collisions can be neglected during the interaction. 
When a* <C a, two-body relaxation can be more efficient 
than Kozai precession in changing stellar orbits (compare 
cqs. 151 and [TU1 and notice that er 2 oc a' 1 at a* <C ro), 
but the number of these stars is negligible. Under these 
conditions, the problem can be tackled by means of re- 
stricted three-body scattering experiments. 

3. SCATTERING EXPERIMENTS 

The integration of the three-body encounter equations 
is performed in a coordinate system centered at the lo- 
cation of M\. Initially, the binary (of mass ratio q and 
eccentricity e) has a randomly-oriented orbit with M2 at 
its pericenter: stars move in the x~ y plane with percen- 
ters along the positive x-axis and random orbital phases. 



The initial conditions of the restricted three-body prob- 
lem problem are then completely defined by 6 variables, 3 
for the binary and 3 for the star: 1) the inclination of the 
orbit of the binary, 9, i.e. the angle between the angular 
momentum of the binary and the z axis; 2) the longitude 
of the secondary hole ascending node, I; 3) the argument 
of the pericenter of the secondary hole, 4> (if e ^ 0); 4) 
the semimajor axis of the stellar orbit, a*; 5) the normal- 
ized (by the angular momentum of a circular orbit with 
the same semimajor axis) angular momentum of the star, 
j*; and 6) the orbital phase of the star, p*. We start each 
scattering experiment by generating 6 random numbers, 
with cos 9 evenly sampled in the range [— 1 , 1] , and both / 
and <fi uniformly distributed in the range [0, 2ir]. We sam- 
ple a* logarithmically around a (the range is described 
in detail below) and j 2 randomly between and 1 (cor- 
responding to an isotropic distribution). Given the of 
a star, we numerical integrate one revolution of a Kcplc- 
rian orbit with eccentricity e* = (1 — J 2 ) 1 ^ 2 , and derive 
p*{t) as a function of time t. Then the initial orbital 
phase for the scattering experiment is drawn from the 
distribution function /(p*) = dt/dp*. 

Having defined the initial conditions, the orbit of each 
star was followed by integrating the coupled first-order 
differential equations 

r = v (12) 
v = - G ±Ml^lA, (13) 

where r and v are the position and velocity vectors of 
the star and rj is the position of the ith (i = 1,2) 
MBH. When e ^ 0, wc included a subroutine to nu- 
merically compute the positions of the two holes at each 
timestep. The units in the numerical computation were 
G = A/12 = a = 1 (with A/12 = A/i + A/ 2 ) and the 
integrator was an explicit Runge -Kutta method of or- 
der 8 (dopri8. lHairer et alJll987[ ). with a fractional er- 
ror per step in position and velocity set to 10~ 13 . The 
integration was stopped if one of the following condi- 
tions was satisfied: (1) the star left the sphere of ra- 
dius a(10 10 g) 1 / 4 , where the quadrupole force from the 
binary is ten orders of magnitude smaller than GA/12/a 2 : , 
with positive energy; the physical integration timescale 
exceeded 10 10 yr; (3) the number of required integra- 
tion timesteps reached 10 8 . Conditions (2) and (3) were 
adopted to save computational time, as a small fraction 
(< 3% depending on q and e) of stars are scattered into 
wide, bound orbits and may survive many revolutions. 
We have tested our cod e by reproducing Figures 4 and 
6 of ISesana et al. (2008) (who used full three-body scat- 
tering experiments), and found excellent agreement. 

4. TESTS 

To understand the dependence of our results on vari- 
ous properties of the MBHB, such as q and e, and the 
impact of general relativistic effects, we have performed 
a number of tests with N = 10 4 stars in each run. The 
initial semimajor axis of the intruder star was sampled 
logarithmically in the interval [l/2a, 2a], where three- 
body interactions are expected to be the strongest. We 
recorded the minimum separation between the stars and 
the holes during each scattering experiment, and ana- 
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Fig. 1. — Close-encounter probability for bound stars interacting 
with a MBHB of mass ratio q = 1/81 and orbital eccentricity 
e = 0.1. The vertical axis shows the fraction of stars N/Ntot with 
closest approach distance r m j n < r. Solid lines: multi-encounter 
cross section. Dashed lines: single encounter cross-section. The 
upper solid and dashed curves refer to the primary hole, the lower 
ones to the secondary. The short vertical lines mark the positions 
of the tidal and Schwarzschild radii of the two holes for M-j = 1 , 
°"loo = 1, a = aoi and 7 = 2. The sharp cutoff at N/Ntot — 10~ 4 is 
an artifact of small number statistics in the scattering experiments. 

lyzed the results in terms of the fraction of stars reaching 
a given distance from a member of the pair. In the case 
of unbound stars, if the initial distribution of pericen- 
ter distances is uniform, such fraction has t he physical 
mean ing of a close-encounter cross section (jChen et al.l 
120081 ). While for the bound stars considered here, the 
concept of cross section no longer strictly applies because 
the initial pericenter-distancc distribution is not uniform, 
for convenience we shall still refer to this fraction as the 
close-encounter cross section in the following. 

4.1. Close-encounter cross section 

We performed scattering experiments for q = 1/81 and 
e = 0.1, where each star was allowed to encounter the bi- 
nary as many times as required before the integration was 
stopped. Then the minimum separation between the star 
and each hole during the entire course of the interaction 
was recorded for the calculation of the "multi-encounter 
cross section" . We also recorded the first minimum sepa- 
ration (a local minimum in the distance-time curve) dur- 
ing the encounter between the star and each hole, and 
calculated the "single-encounter cross section" . The lat- 
ter can be viewed as the interaction probability in the 
case of an isolated MBH. The resulting cross sections are 
plotte d in Figure [1] As already shown by iChen et al.l 
<|2009f ) . the multi-encounter cross section for Mi is dra- 
matically higher than the single-encounter probability: 
when r ~ 10~ 4 a, corresponding to the tidal radius of a 
primary hole with M7 = 1 embedded in an isothermal 
cusp with (J100 = 1, at separation a = ao, the multi- 
encounter cross section of M± is nearly 3 dex larger than 
that for single encounters. In the following, we shall refer 
to the multi-encounter cross section at r t \ as the "tidal 
disruption cross section". Because this is much larger 
than the corresponding cross section for the secondary, 
stellar disruptions by an unequal-mass binary will be 
dominated by the primary hole, and will be the focus 
of our analysis. 

4.2. Dependence on q and e 
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Fig. 2. — Multi-encounter cross section for bound stars interact- 
ing with the primary hole of a MBHB with e = 0.1 and different 
mass ratios q. The short vertical lines mark the locations of the 
tidal radii of M\ for the same parameters used in Fig. [T] 
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Fig. 3. — Same as Fig. [2] but for q = 1/81 and different values 
of the binary eccentricity e. 

To study the dependence of the close-encounter cross 
section on the binary mass ratio, we performed three 
additional sets of scattering experiments for e = 0.1 and 
9 = 1/9, 1/243, and 1/729, each using 10 4 stars. The 
results show (Figure [2J that, as q increases, the multi- 
encounter probability decreases from 41% (q = 1/729) to 
1.1% (q = 1/9). This is because, as the perturbing force 
from the secondary hole becomes stronger, a star is more 
susceptible to ejection. Figure [3] shows the dependence of 
the multi-encounter cross section on binary eccentricity 
at fixed q = 1/81. For r/a > 10~ 6 , the cross section 
varies at most by a factor of 3 as e increases from 0.1 to 
0.9. When r t \ja ~ 10~ 4 , the tidal disruption probability 
does not depend significantly on eccentricity. 

5. PROPERTIES OF DISRUPTED STARS 

To understand the physical processes responsible for 
the enhancement of the tidal disruption probability, we 
need to investigate the properties of the disrupted stars. 
We performed new scattering experiments aimed at cov- 
ering the whole parameter space of the interacting stel- 
lar population, extending the range of semimajor axis 
a* from [a/2, 2a] to [a/20, 20a]. We ran four sets of nu- 
merical experiments, each consisting of 5 x 10 4 stars, for 
varying binary eccentricities and q = 1/81. A star is 
counted as disrupted if its separation from the primary 
hole becomes smaller than r t \. (To calculate r t i/ao, the 
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fiducial parameters M-j 
used.) 
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5.1. Phase-space distribution 

The semimajor axis of a MBHB typically shrinks by 
a factor of ~ 10 during the process of cusp erosion via 
three-body scatterings (jSesana et al.1 120081 ) . Below we 
scale the same scattering experiments and present re- 
sults for two cases, a = ao and a = do/ 10. Figure [4] 
shows the distribution of disrupted stars in the a» — j~ 
plane for e = 0.1. The fraction of disrupted stars ex- 
ceeds 19% in the case a = a /10, and is close to 13% 
for a = ao- Many of the stars that get disrupted are 
initially located outside the tidal loss cone, showing that 
j* is not conserved during the three-body interaction; on 
the other hand, stars initially outside the interaction loss 
cones get disrupted only rarely. Figure 0] also shows an 
excess of stars at the resonance radii a* = a(m/n) 2 / 3 , 
where m, n = 1, 2, 3..., indicating the importance of res- 
onant interactions in refilling the tidal loss cone. 

For a better understanding of the nature of disrupted 
stars we depict in Figure [5] their distribution in the 
a* — jz* plane. Both theoretical and numerical stud- 
ies show that for stars that lie inside the binary orbit 
(with semimajor axis a* < a/2), the angular momentum 
component parallel to the binary orbital angular momen- 
tum, j z *, does not change, while the angular momentum 
comp onent perpendicular to j x * undergoes secular evolu- 
tion (jKozai lfl962t iLidov II1962I ). This implies that stars 
in the wedge- like region \j z *\ < j\ c {rti I a*) will undergo 
secular evolution and finally enter the tidal loss cone and 
get disrupted. Figure confirms that the majority of the 
disrupted stars with a* < a/2 have \j z *\ ^ iic(?"ti/a*), 
i.e. lie within the region delimited by the solid lines. 
When a* > a/2, however, stars on eccentric orbits cross 
the orbit of the secondary hole, and can get disrupted 
even if \j x +\ ^> j\ c {rti / a*) ■ These stars are difficult to 
model as their orbits are chaotic. The size of the interac- 
tion loss cone relative to the Kozai wedge increases with 
a (r t i/a decreases). As a result, strong chaotic three- 
body interactions rather than cumulative secular effects 
are responsible for the majority of the disruptions. 

For an isotopic stellar distribution, the fraction of stars 
having semimajor axis in the range (a*, a* + Aa*) that 
reside inside the "Kozai wedge" is given by 

rjlc[(l+e)a/a,] r j lc (r t i /a, ) 

fx {a* /a) = J dj» / dj z * 

= 2ji c (r t i/a*)j lc [(l + e)a/a»] - 2jf c (r tl /a t ). (14) 

In a binary system with Mf = 1, e = 0.1, 
cioo = 1, 7 = 2, and a = ao, the mean frac- 
tions /if(a*/a) in the strong-interaction regime a/2 < 
a* < 2a are (0.0089,0.025,0.044,0.074) for q = 
(1/9,1/81,1/243,1/729). These theoretical estimates 
are significantly smaller than the tidal disruption proba- 
bilities derived in § 14.21 except when q = 1/9, highlight- 
ing the importance of chaotic interactions in the rcpop- 
ulation of the loss cone. For q = 1/9, the theoretical 
tidal disruption cross section becomes comparable to the 
numerical one, because stars in the chaotic-interaction 
regime are more susceptible to early ejection. 



Figure [5] shows the distribution of disrupted stars in 
the a* — j z » plane for the extreme e = 0.9 case. Note 
that the interaction loss cone is Jz c [(l + e)a/a»] when e 
is large. The number of disrupted stars increases sig- 
nificantly relative to the low binary eccentricity case, to 
about 38% when a = ao/10 and 24% when a = ao. This 
enhancement occurs as more stars now cross the orbit of 
the secondary hole and interact chaotically with the bi- 
nary. Moreover, any trace of the Kozai wedge disappears 
in this case. This is because, when e = 0.9, the apoastron 
of the secondary hole is 0.1a; therefore, even stars with 
a* <C a experience strong interactions with the secondary 
hole that destroy the secular, coherent accumulation of 
the Kozai mechanism. 

5.2. Disruption timescales 

During a scattering experiment a bound star may en- 
ter the tidal sphere of the primary black hole many times 
before it is ejected. When calculating the disruption 
rate, it is the time when the star first crosses r t \ that 
is relevant: in the following, we refer to this time as the 
"tidal disruption timescale" . In our numerical integra- 
tions, we record the times when each star first reaches 21 
different separations logarithmically distributed within 
the range [log(100r t i/ao), log(r t i/100ao)]. Then, for any 
binary separation between ao/100 and 100ao, the tidal 
disruption timescale can be derived by interpolating be- 
tween the recorded times. For example, the time when 
a star reaches 10r t i/ao can be viewed as the tidal dis- 
ruption timescale (in units of the binary orbital period) 
when the binary has shrunk to a = ao/10. Figure [7] 
shows the tidal disruption timescale, r (in unit of the 
binary period P), as a function of the initial semimajor 
axis, for e = 0.1. The dashed line indicates the stel- 
lar orbital period, P, = P(a»/a) 3 / 2 . The cross symbols 
clustered around the dashed line represent stars that are 
disrupted within one orbital period because their initial 
pericenter distances are smaller than the tidal radius of 
Mi. The solid lines in Figure [7] trace instead the Kozai 
timescale. The formula derived in equation (|10|) applies 
to stars that orbit close to the primary hole while the sec- 
ondary is far away ("inner problem"). In this case, the 
quadrupolc force exerted by the secondary on the star, 



Fn 



(Gmtl/ 2 a»)a , 



(15) 



causes a torque a* Ft that changes the angular momen- 
tum of the star on the timescale 



i J* 
K ~ dJ*/dt 



m^GJVha*) 1 / 2 
a*Fr 



P 

2nq \a„ 



3/2 



(16) 



In spite of the many simplifications in the derivation of 
equation (fT6"|) . T' K differs from the actual Tk by only a 
factor of 4/3. 

Based on our understanding of the Kozai effect in the 
inner problem, we can now estimate the Kozai timescale 
for the case of orbit crossing between the intruder star 
and the secondary hole ("outer problem"). When a* 3> 
a, the quadrupolc force exerted on the star by the binary 
is 

and a*i*T is the corresponding torque. Since the peri- 
center of a star must be smaller than a (star lies in the 
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Fig. 4. — Distribution of disrupted stars in the a* — j% plane, assuming e = 0.1. Left panel: a = ao/10. Right panel: a = ao- The solid 
and dashed lines delineate, respectively, the tidal loss cones and the interaction loss cones. The parameters of the binary-stellar cusp are 
q = 1/81, My = 1, o"ioo = 1, and 7 = 2. The excess of disrupted stars at radii a„ = a(m/n) 2 / 3 , where m, n = 1, 2, 3..., is due to resonant 
interactions. 
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Fig. 5. — Same as Fig.[4]but in the a* — j z * plane. 
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Fig. 6.— Same as Fig.[5]but for e = 0.9. 
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interaction loss cone) for it to be scattered into the tidal 
loss cone (see § 15.11) , the maximum stellar angular mo- 
mentum is m*ji c (a/a*)y/ GM% a* . The Kozai timescale 
in the outer problem is then given by the ratio of the 
maximum angular momentum and the torque, 

m„ji c (a/a*)(GM 1 a Ji ) 1 / 2 2 
T K oc tx (a*/a) P. (18) 

Since the transition radius between the inner and outer 
problems is a* ~ a/2, continuity between equations (pi 
and (HI]) yields 



T K = 




(a* < a/2) 
(a* > a/2) 



(19) 



In Figure [7] the crosses mark stars initially inside the 
Kozai wedge, so that crosses around the solid lines rep- 
resent stars fed into the tidal loss cones via the Kozai 
mechanism. The dots mark instead stars initially out- 
side the Kozai wedge; these participate to the chaotic 
loss-cone repopulation, and are disrupted on timescales 
that are typically longer than the Kozai timescale. As 
the binary orbital separation increases, less stars are dis- 
rupted by the Kozai mechanism and stars on chaotic or- 
bits make a larger contribution to the repopulation of the 
tidal loss cones. 

Figure [5] shows the tidal disruption timescales for the 
high-eccentricity, e = 0.9 case. The cross symbols no 
longer cluster along the solid lines and the dispersion 
is larger; this is because resonant interactions are now 
stronger and partially suppress the secular Kozai evolu- 
tion. 

5.3. Relativistic and cusp-induced apsidal precession 

Relativistic effects close to the central MBH and the 
presence of an extended stellar cusp can both induce 
the precession of the pericenter of a stellar orbit. These 
would act to suppress the Kozai mechanism if the asso- 
ciated precession rates were higher than the Kozai rate. 

The effect of relativistic precession can be understood 
as follows. The pericenter of a star with semimajor axis 
a* and eccentricity e* undergoing a secular Kozai evolu- 
tion precesses at an approximate rate 



U1 K 



(l-e2)V 



Tn TK ^ a *y 



(20) 



When e* grows to e* ~ 1 — rti/a*, the star gets tidally 
disrupted by the primary hole during the pericenter pas- 
sage. Using (1 — e*)a* ~ 2r t \ and equation (|T9|) . equa- 
tion (|2T)|) gives the Kozai precession rate at the tidal dis- 
ruption radius, 



. i7§^ (^r 1/2 (t) 2 ( a *< a / 2 ) 
M> (^r 1/2 (tr 3/2 («.>«/2) 



.(21) 



The general relativistic precession rate is instead 

Q-kGMi , 37r / rsx\ /a*\ -3 / 2 



WGR 



(1 - el)c 2 a, 



2P(a) \r t 



The condition wqr = admits only one solution, 
a*,cri = (8v / 2£) _2/ ' 7 a, where 
1/2 

(23) 



16 \ r si/ \ r *i 



~1.64x itf^-ifftM^o^qV-WQ-^ W(24) 

If a* < a* lCr ii the Kozai evolution is suppressed by GR 
precession. Note that when £ < 1, wgr > ojk, and 
stellar disruptions from the secular Kozai mechanisms 
are suppressed for the entire stellar population, leaving 
only chaotic encounters to contribute to the tidal dis- 
ruption rate in this regime. Since the ratio between the 
Schwarzschild radius and the tidal radius of an MBH, 

r s /r t ~ 0.19M 7 2/3 (R o /r*)(M,/M Q ) 1 / 3 , (25) 

increases with hole mass, GR effects typically become 
important when Mbh > 3.6 x 10 6 M Q (i.e. r t < 10rs for 
solar type stars). 

Figure [5] shows the quantity a* iCr i/a as a func- 
tion of q for different combinati ons of the parameter s 
(7, a /ao, M 7 ), assuming M 7 cx cr 4 (jTremaine et al.ll2002j) . 
Following equations (|2~H) and ((9J , the curves move toward 
the upper-right direction as r t i/rgi decreases or r t \/a in- 
creases. In the upper-right corner of each curve in the 
q — a*,cri/a plane, the Kozai mechanism is effective, while 
in the lower-left corner wgr > ^>K and stellar disruptions 
are suppressed. When q > 0.01, the GR precession does 
not significantly affect the Kozai evolution of stars in 
the strong-interaction regime (a* £ [a/2, 2a]), when the 
stellar-disruption fraction is the highest (see Fig. ITO]) . 
F° r <7 ^5 0.01, however, GR effects are important in the 
strong-interaction regime, especially in the case of steep 
stellar cusps (7 > 2), compact MBHBs (a/ao < 0.3), or 
massive primary holes (M7 > 3). 

To simulate numerically the effects of GR, we have 
run a series of scattering experiments for different val- 
ues of q and e using the p seudo-Newtonian potential of 
iPaczvhski fc Wiit"aT(|1980l ) . We integrated the equations 



of motion 



r = v 



M»(r-ri) 



(26) 
(27) 



(22) 



where r$i is the Schwarzschild radius of the ith hole. 
Each set of experiments followed 10 4 particles sampled 
in the range a* E [a/20, 20a]. For illustrative purposes, 
we calculated rsi/a assuming a = ao, M7 = 1, crioo = 1, 
and 7 = 2. For this parameter choice, rn/rsi — 5. The 
integration of the stellar orbit was stopped at l.Olrsj to 
avoid the singularity at the Schwarzschild radius. Fig- 
ure |TD] compares the fractions of the disrupted parti- 
cles in the GR versus the non-GR simulations. When 
q = 1/9,1/81, the suppression of stellar disruptions by 
relativistic precession is important for a* < a* iCr i. When 
q = 1/243, wgr is always larger than cok (£ = 0.43), and 
tidal disruptions are suppressed over the entire range of 
a*. 

The impact of an extended stellar cusp on the Kozai 
evolution can be addressed using similar arguments. 



8 



Chen et al. 




Fig. 7. — Tidal disruption timescales in unit of the binary period P for a = ao/10 (left panel) and a = ao (right panel), as a function of 
the initial scmimajor axis of the intrud er st ars. The dashed line represent the initial stellar orbital periods, and the solid line marks the 
analytical Kozai timescale given by eq. {X9j . The red crosses refer to disrupted stars initially inside the Kozai wedge, \jz*\ < Jic( r *ti/ a *)- 
System parameters as in Fig. [4] 




Fig. 8. — Same as Fig.[7]but for e = 0.9. 
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Fig. 9. — Tracks of a* jC ri/fi as a function of q for different values of 
7 (top), a/ao (middle), and M7 (bottom). The variables are labeled 
in the top-right corner of each panel, and the fixed parameters in 
the lower-left corner. 
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The cusp- induced precession rate is (jlvanov et ID [2005; 
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are dominant. 

6. DISRUPTION RATES FOR STATIONARY BINARIES 

Having set the properties of the disrupted stars and 
tested the limitations of our approximations, we can now 
proceed to the calculation of the disruption rates. We 
fix the orbital separation a and the eccentricity e of the 
MBHB, and assume that the stellar cusp surrounding 
Mi is isotropic and composed of solar-type stars only. 
The initial stellar distribution function, fo(a>*,j*,jz*) = 
duo/da*, can then be written as 



Fig. 11. — Tracks of j% cri as functions of a* for different com- 
binations of the parameters 7 and a/ag. Stellar orbits satisfying 
lu c < Cjk lie below the tracks. 

where K = (0.5, V^/n) for 7 = (2, 1.5) respectively, and 
M*(a*) is the stellar mass enclosed within a sphere of 
radius a*. The dependence of tb c on (1 — e 2 ) 1 / 2 indi- 
cates that the stellar cusp affects mostly circular orbits. 
Assuming the broken power-law distribution described 
in § [3J the condition io c = lok corresponds to a critical 
angular momentum 



7 ' 2 .: 

J *,cn 



(ls/ie)/^ 1 



7-3 



(if) 

7-3 



V a 



above which the Kozai effect is suppressed by cusp- 
induced precession. Figure [TT] shows tracks of j 2 cri in 
the a, — jl plane: once again, the phase space where 
the Kozai mechanism can operate is greatly reduced by 
cusp-induced precession. The suppression is expected to 
be more significant if a or 7 increase, as the stellar mass 
enclosed by the stellar orbits increases. 

We mimicked the effect of a broken power-law stel- 
lar cusp (§ [2]) by including an external potential in the 
equations of motion. We set 7 = 2 to maximize the ef- 
fect of the cusp (see Fig. [TT]), and ran two sets of 10 4 
non-GR scattering experiments for e = 0.1 and 0.9 re- 
spectively. The binary parameters were set to q = 1/81, 
e = 0.1, M7 = 1, trioo = 1, and a = ao- Figure [T2l 
show the initial values of a*,_7 2 for the disrupted stars, 
while the fraction of tidal disruptions as a functions of 
a* is shown in Figure [T51 When e = 0.1, tidal disrup- 
tion is almost completely suppressed when a* < a/2 for 
stars with j» > j* tCT i] when a* > a/2 a small fraction 
of stars with j* > j*, C ri still get disrupted from chaotic 
interactions. When e = 0.9, the suppression of stellar 
disruptions is also appreciable: a large fraction of stars 
with j* > j* ;C ri are disrupted, however, because their 
orbits undergo chaotic intersections with the highly ec- 
centric orbit of the secondary MBH. The dotted lines in 
Figure [TBI show that, when e is small, the tidal disrup- 
tion fraction in the presence of a cusp potential is ap- 
proximately equal to the fraction of disrupted stars with 
j* < j*,cri when the cusp is neglected, with an error that 
is less than a factor of two. When e is large, however, this 
method severely underestimates the true disruption frac- 
tion because chaotic three-body interactions at j* > j* )Cr i 



da* 



2(3-7) M 2 



"0 



Mr. 



2-7 



(30) 



where dno /da* is the number of stars per unit semima- 
jor axis a*, and the right hand side follows from the 
definition of ao- We do not consider any sharp cutoff 
in the stellar distribution caused by stellar collision at 
small radii. Collisions will likely result in a shallower 
inner density profile rather then a well defined cutoff 
(jFreitag fc Bend [2002). Moreover, the inner cusp may 
be repopulated by efficient g as inflow-in duced star forma- 
tion during galaxy mergers (Zier 2006]). Here we assume 
7 = 2 and 1.5 to account for this uncertainty. 

The stellar disruption rate at time t can then be cal- 
culated from 



(a* < a/2) 

(29,) 
(a* > a/2) 



F(a*/a,t)(dn /da*)da*, 



(31) 



where F(a*/a,t)dt is the fraction of stars with semi- 
major axis a* that are disrupted in the time interval 
(t, t + dt). If the loss-cone refilling is entirely due to the 
Ko zai effect, then F(a */a,t) can be derived analytically 
as (|Ivanov et al.ll2005f ) 



F(a*/a,t) 



fK(a*/a) 
T K (a*/a) 



exp(-t/T K ), 



(32) 



where fx and Tk are given by equations (TTT)) and (TTJ 
The chaotic nature of strong three-body scattering events 
prevent the possibility of a simple analytical description; 
therefore, the total disruption rate, including the con- 
tribution of resonant interactions, has to be computed 
numerically from scattering experiments. We divide the 
(a*/20a, 20a*/a) interval into 52 equal logarithmic bins. 
Given the parameters e and r t \/a, we derive the function 
Fi(r) numerically for each a* /a bin using the recorded 
tidal disruption timescales in the corresponding scatter- 
ing experiments, so that i 7 i(r)Ar is the fraction of stars 
in the ith bin that are disrupted in the time interval 
(r, r + At). When deriving Fj(r), stars sampled in the 
range j*(ai) < ji c {rn/ai) are excluded. The total stellar 
disruption rate at time t is then given by 



52 



F i (t/P)(dn /da*){a i )Aa l 



(33) 



i=i 



where and Aa^ are the central semimajor axis and the 
width of the ith bin. To calculate the numerical stellar 
disruption rate, one must specify P and dno/ da* in phys- 
ical units. From the definition of ao and equations pip 



10 



Chen ct al. 





Fig. 12. — Distribution of disrupted stars in the a, — jl plane when the potential of a stellar cusp is included (see text for the parameters 
of the scattering experiments). The eccentricity of the binary is set to e = 0.1 {left) and 0.9 {right). The solid lines show the tracks of j% .. 
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Fig. 13. — Fraction of disrupted stars as a function of a* in scattering experiments with {solid curves) and without {dashed curves) 
the cusp potential. The dotted lines show the disrupted fractions estimated according to the j t < cr ; criterion. The left (right) panel 
assumes e = 0.1 (e = 0.9). 



and (1301). we have 



^ M 7 , ..^3/2 ( ± 

a 



P ~ (3 x 10 b yr) Q (3 - 7 ) ; 



'100 



3/2 



with Q = q^^-^/il + q) 1 / 2 , and 



2-7 



(34) 



(35) 



with a = (2 — 7)/(3 — 7). Figure [14] shows the total 
stellar disruption rates calculated from equation (|33j) for 
a MBHB with a = ao and a = ao/10 (and the standard 
system parameters q = 1/81, e = 0.1, M 7 = 1, crioo = 1, 
and 7 = 2). The stellar disruption rate remains constant 
for a timescale P/q before decreasing rapidly, consistent 
with the Kozai time scaling. Compared to the rates 
for single MBHs fed by two-body relaxation, typically 
10~ 4 — 10~ 5 yr -1 , the rates on the plateau are orders 
of magnitude higher. As a decreases, the plateau rate 
increases while its duration shortens. Since the Kozai 
timescale increases as a 3 / 2 and the number of stars en- 
closed in the Kozai wedge scales as a 1 / 2 , if a shrinks by 




10" 10" 10" 

t [yr] 

Fig. 14. — Stellar disruption rates as a function of time for a 
stationary MBHB with semimajor axis a = ao/10 {solid curves) 
and a = ao {dashed curves). The thick lines are the numerical 
rates derived from scattering experiments, and the thin lines are 
the analytical rates. System parameters are the same as in Fig. [4] 
and yield P ~ (400, 13) yr for a = (ao, ao/10). Fluctuations in the 
numerical rates at early times are due to poor statistics. 
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FlG. 15. — Numerical (thick curves) and analytical (thin curves) 
stellar disruption rates as a function of time for different combina- 
tions of the parameters e and 7 (see text for details) . The scmima- 
jor axis of the binary is a = ao, and all other system parameters 
are the same as in Fig. [4] 

a factor of 10 the plateau phase becomes a factor 10 3 / 2 
shorter, while the disruption rate at the plateau increases 
by a factor of 10. The figure also depicts the analytical 
disruption rates calculated with equation (|31[) for com- 
parison. These agree very well with the numerical results 
during the plateau phase, indicating that the tidal loss- 
cone refilling is initially dominated by the Kozai mech- 
anism. At later times, stars inside the Kozai wedge are 
mostly depleted, and the analytical and numerical rates 
start to deviate from each other. Deviations in the post- 
plateau phase increase with increasing binary orbital sep- 
arations. According to our numerical calculations, about 
(1.8 x 10 4 , 1.1 x 10 5 ) stars with a/20 < a, < 20a are dis- 
rupted over 10 s years by binaries with a = (ao/10, do). 
The corresponding numbers in the analytical approxi- 
mation are (1.1 x 10 4 ,3.5 x 10 4 ). The difference high- 
lights the importance of close, resonant encounters with 
the secondary hole, which change the stellar orbits in a 
chaotic manner and fuel the tidal loss cone. Figure [T5] 
shows the dependence of the disruption rate on the pa- 
rameters e and 7. Increasing the binary eccentricity only 
affects the rate in the post-plateau, chaotic-interaction- 
dominated phase. 

The above results show that, whereas chaotic scatter- 
ings dominate the total number of disrupted stars, the 
Kozai theory provides a reasonably good description of 
the disruption rate during the plateau phase, as well as 
the correct order of magnitude of the total number of dis- 
ruption. We can then use the Kozai theory to predict the 
scaling of the plateau rate with the system parameters, 

,v M*f K ,„„s 
~T~' ' ' 

Here M* is total mass of the interacting stars, M* oc 
Mig(a/ao) 3 ~ 7 (from the definition of ao), fx is the frac- 
tion of stars in the Kozai wedge, fx oc (rti/a) 1 / 2 , and 

Tk is the Kozai timescale, TV oc q~ 1 a 3 l 2 M x 1 ' 2 (assum- 
ing that a* ~ a). Substituting into equation ()36|) . and 
using the definitions of r tl and ao, we finally obtain in 
the limit q <C 1: 

JV,oc(3- 7 rV 4 - 2 ^ A3 -^ 1 M- 1/3 at 



Fig. 16. — Stellar disruption rates as a function of time for a 
stationary MBHB with semim ajor axis a = ao. Binary parameters 
are the same as those in Fig. 1141 The solid line is derived from 
the scattering experiments where both GR and cusp effects are 
included, while the dashed line from the experiments when both 
are neglected. The dotted line is also derived from no-GR/no-cusp 
experiments, but only using stars with a* > a» cr i and jr* < j*, cr i. 

oc (3 - 7 )~ Y 4 ^)/(3-7) f^j M V3 ) (37) 

wher e we have adopted Mi oc l|Tremaine et al.l 
l2002f ) in the second proportionality. According to equa- 
tion (|3"T)) . when q = 1/81, as 7 decreases from 2 to 1.5, 
the peak stellar disruption rate drops by a factor of 40, 
consistent with the rates in Figure Hoi If 7 = 2, the dis- 
ruption rate is proportional to a^ 1 , consistent with the 
left panel of Figure Q31 The above analysis also implies 
that, when the stellar density profile is as steep as 7 ~ 2, 
the peak stellar disruption rate is not sensitive to q, as 
shown by IChen et al.l (|2009f ) . Assuming the M\ — cr» re- 

lation, the peak rate should be proportional to M l ' . 

To investigate numerically the impact of GR and cusp- 
induced precession on the stellar disruption rate, we ran 
an additional set of 10 4 scattering experiments that in- 
cluded the two effects simultaneously (sec § 15.31 for de- 
tails). The binary parameters were set to q = 1/81, 
e = 0.1, 7 = 2, My = 1, ctioo = 1, and a = ao- The 
resulting stellar-disruption rate is shown in Figure [T|)] 
with the solid line, and is compared to the case in which 
the two effects are neglected (the dashed line). Because 
chaotic scatterings are not suppressed by secular effects, 
the two curves differ by only about a factor of two. The 
contribution from chaotic scatterings in the GR+cusp 
experiments can be gauged by the difference between 
the solid and the dotted curves, the latter derived by 
considering only stars with a* > a* jCr j and < j* )Cr i 
in the no-GR/no-cusp experiments. For binaries with 
larger g, larger e, smaller 7, or larger a/ao, stars in the 
chaotic-interaction regime contribute more to the dis- 
rupted stellar fraction, and the GR/cusp- induced sup- 
pression of tidal disruptions is milder. Disruptions rates 
derived by scattering experiments that neglect GR and 
cusp effects, together with the analytic scalings given by 
equation (|37|) , can therefore be considered valid for bina- 
ries with q > 0.01. 

7. DISRUPTION RATES FOR DECAYING BINARIES 
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Recent calculations based on scattering experiments 
and ignoring stellar disruptions have suggested that 
both the binary scmimajor axis and eccentricity will 
evolve rapidly du ring three-body int eractions with ambi- 
ent bound stars (jSesana et al.l 120 08). On the one hand, 
in a shrinking MBHB the interaction loss cone, the tidal 
disruption timescale, and the cusp stellar distribution all 
change with time, and this evolution will affect the stellar 
disruption rate. On the other hand, tidal disruptions halt 
the exchange of energy and angular momentum between 
the stars and the binary, altering its dynamical evolution. 
In order to compute a more realistic tidal consumption 
rate, a hybrid model is required that takes into account 
stellar ejections as well as disruptions, and that solves 
for the time evolution of the binary-stellar cusp system. 

In evolving MBHB systems, the ratios rsi/a and 
M*(a)/Mi vary with time. In this case, simulating GR 
and cusp effects becomes extremely time consuming, be- 
cause additional scattering experiments need to be car- 
ried out whenever rgi/a or M*(a)/Mi changes. For this 
reason, in this section we do not consider GR and cusp 
effects, and use the Newtonian scattering experiments 
without cusp potential to calculate the stellar disruption 
rate. We restrict the following discussion to the case 
q > 1/0.01, where GR and stellar cusp precession sup- 
press the stellar disruption rate by only a factor of two. 
When q <C 0.01, our test calculations with q = 1/729 
show that GR and cusp effects can suppress the stellar 
disruption rate by as much as a factor of ten. 

7.1. Fate of interacting stars 

Figure [17] compares the fraction of stars that arc 
ejected from the system with those that are disrupted 
in scattering experiments with q = 1/81. Tidal disrup- 
tions produce two interesting effects: (1) when a» ~ a a 
significant fraction of stars experience strong interactions 
with the secondary hole and cross the tidal radius of M\ 
before being ejected; (ii) frequent tidal disruptions occur 
even when a* <C a, a regime where ejections arc rare. 
These are stars that are driven into the tidal loss cone 
by the Kozai mechanism. Tidal disruptions have then 
the double effect of partially suppressing stellar ejections 
(especially when the binary eccentricity is large) and at 
the same time of extending inward the influence domain 
of the binary (the a* /a interval where the black hole pair 
can alter the stellar cusp). Figure fT8l shows the distribu- 
tions of changes in specific energy and angular momen- 
tum (z-componcnt) for the stars that are ejected and for 
those that are disrupted. Such distributions are narrowly 
peaked around zero in the case of the disrupted popula- 
tion, while are much broader and skewed towards positive 
values for the ejected component. The evolution of the 
MBHB is then determined by stellar ejections, since on 
average the disrupted stars do not exchange energy and 
orbital angular momentum with the binary. 

7.2. Hybrid model 

We finally describe our hybrid model. Given the 
binary-cusp system parameters q, 7, M7, and ctioo, we 
calculate the corresponding value of the initial binary 
scmimajor axis oq. We use a* to describe the absolute 
scmimajor axis of interacting stars, and only consider the 
relevant portion of the cusp enclosed in the a* interval 



[10 _3 a , 100a ]. This range is binned in 100 equally log- 
spaced bins labelled by the index i, and the initial stellar 
mass in each bin is given by m, = m^Aa* idno(a* i)/da* 
(i = 1,2, 3,..., 100), where the ccntroid of the ith 

bin and Aa*,, is the bin width. At t = the MBHB is at 
separation ao with eccentricity eo and period Pq ; we then 
evolve the system numerically forward in time according 
to the equations 



ejfe+i =efe 



AE k 
1 - 



-Oft, 



( AE k , 2AJ A 



2e fc V E k 



(38) 
(39) 



where the index k (fc > 0) labels the timestep, a k and e k 
are the binary semimajor axis and eccentricity, Ek and Jk 
are the energy and angular momentum of the binary, and 
A refers to the variation in the fc-th timestep Atk- The 
increments AJk and AEk, depend on the mass ^ Am^fc 
that interacts with the binary in the k-th timestep. The 
subtlety lies in properly extracting Am^j. from the set of 
scattering experiments described in Section 3. Numerical 
experiments are carried at a fixed binary separation, and 
the relevant parameter in determining the fate of a star 
is the ratio s = a* /a. In our runs we sample the range 
1/20 < s < 20, and this interval is divided in equally 
log-spaced bins labelled by the index j as Sj. For each 
Sj bin, we construct the functions df/dr\j, d£/dr\j and 
dj/dr\j, which are the differential fractions of ejected 
stars, mean energy exchange, and mean angular momen- 
tum exchange as a function of r, the time expressed in 
unit of the binary period. The trick is to assign to each 
bin a*^ the right Sj value as the binary semimajor axis 
a evolves, and to properly connect the physical time t 
describing the evolution of the system to the 'scattering 
experiment time' r (expressed in units of P). For the 
time being, let us ignore, for simplicity, the eccentric- 
ity evolution. The integration scheme then proceeds as 
follows. 

Consider the first timestep At®. In each a*^ bin, the 
amount of ejected (or disrupted) masa3 in this first 
timestep is 



Ami 



(It 



Jo 



^0 



(40) 



where jo identifies the Sj bin to which the a*,, stellar 
bin belongs in the first timestep. If a particular a*^ bin 
lies outside the 1/20 < s < 20 range, then it does not 
contribute in the evolution of the binary at the consid- 
ered timestep. The binary separation a is then evolved 
according to equation (|3"5)) . where 



AE 



E 



d£_ 

rf7 



(r = 0) 



At 



Ami 



(41) 



is given by summing over all a*,,. We accordingly shrink 
the binary from ao to a\ . 

Consider now the second timestep At±. Since the 
stellar distribution changes with time, in principle one 

6 Here we do not distinguish between ejected and disrupted stars. 
The distribution df/dr is actually split in df/dr e j and df/dT^ a to 
account for both components. 
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Fig. 17. — Fractions of disrupted (solid lines) and ejected (dashed lines) stars as a function of stellar semimajor axis a*, for e = 0.1 
(left panel) and 0.9 (right panel). Black thin lines: a = ao/10. Red thick lines: a = ao. The dotted line shows the ejection fractions if 
disruptions are not taken into account. All other system parameters are as in Fig. [4] 




Fig. 18. — Normalized distributions of changes in specific energy 
momentum (right panel) for the ejected (dashed curve) and disrupted 
angular momentum in unit of (GM120) 1 ' 2 . Line styles as in Fig. 1171 

should carry out new scattering experiments according 
to the updated stellar distribution to derive df /dr\j at 
every timestcp. However, as long as the stars depleted 
during the previous steps are appropriately excluded, the 
original scattering experiments can still be used. For the 
stars in a Sj bin, the elapsed scattering-experiment time 
Tj.i during the first timestep can be solved from the im- 
plicit equation 



- 1 df_ 

dr 



dr = Ami 0) 



(42) 



where df/dr\j is the same function as in the first 
timestep, and jx identifies the new Sj bin to which the a*,, 
stellar bin belongs in the second timestep. In the second 
timestep, the time zero point of the function df/dr\j 1 
shifts to t = Tj t i to exclude the stars with depletion 
timescales shorter than Tj t x, so the interaction mass be- 
comes 



Am,, 1 



4f 
dr 



At\ ( a\ 



.11 



-3/2 



, (43) 



(left panel) and in the z-component of the specific orbital angular 
(solid curve) stars when e = 0.1. Energy is given in unit of GA/12/a, 



where (ai/ao) -3 / 2 accounts for the change in the period 
of the binary as it shrinks from ao to a\. We then evolve 
again the binary according to equation (|38[) . where now 
AE\ is given by 



AE 1 = 



d£ 



Ah ( a x 



31 



"0 \ao 



-3/2' 



For a generic timestep Atk , the interacting mass is 

-3/2" 



Am^i. 
(44) 



Ami k 



dr 



j 1, 



( T = T 3,k)-^- — 



(45) 



where j k identifies the Sj bin to which the a,^ stellar bin 
belongs in the fc-th timestep, and Tj t k labels the value of 
t that solves the implicit equation: 



dr 



dr 



3k 



fc-1 

E 

k'=0 



Ami 



k' ■ 



(46) 



The binary is then evolved according to equation (|38[) . 
where AE k is given by replacing the index 1 with k in 
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Fig. 19. — Evolution of binary separation a (top panels) and 
eccentricity e (bottom panels) for a MBHB with different combi- 
nations of the initial e, 7 values. The curves in the left panels 
are computed neglecting tidal disruptions, while those in the right 
panels include the effect of stellar disruptions. The other system 
parameters are q = 1/81, M7 = 1, and crux) = 1- 

equation (|44|) . 

In this way we account for the fact that, in each stel- 
lar bin a*^, the interacting fraction in the timestep k is 
governed by the stars left in the bin at that timestep fol- 
lowing the ejection occurred in the previous timesteps; 
and that the ejection occurs at a rate given by the Sjk 
bin to which the a*^ stellar bin belongs at the fc-th 
timestep. When we also consider the binary eccentric- 
ity evolution, we interpolate the values of df/dr\j and 
d£ / dr\j between the different eccentricities sampled by 
the scattering experiments e = (0.1,0.3,0.6,0.9); the ec- 
centricity is evolved according to equation (|39j). where 
AJk is computed from the analogous of equation (|44|) . 
but using dj /dr. We only considered the variation in 
the z-component of the angular momentum, because in a 
spherical stellar cluster t he rotational B rownian motion 
of a MBHB is negligible (|Merritt 1 120021 ). 

Our hybrid approach relies on an adiabatic approxi- 
mation, i.e. the MBHB orbital evolution is assumed to 
be slower compared to the typical star-binary interaction 
timescale. This is certainly true for chaotic encounters, 
but it is not so for secularly evolving stars. To justify 
our evolution scheme, we have therefore run test scat- 
tering experiments in which the MBHB was evolved by 
hand, according to the shrinking rates derived with the 
hybrid model. The initial conditions of the test experi- 
ments were the same as in § [3] and the stellar-disruption 
rates were calculated following the scheme described in 
§[6] This setup allowed us to directly measure the disrup- 
tion rate caused by an inspiralling binary on a population 
of stars drawn from a chosen density distribution. Fig- 
ure [20] compares the stellar-disruption rates derived in 
this fashion to those given yielded by the hybrid model; 
the agreement between the two is quite good, validating 
our orbital integration scheme. 

7.3. Results 

Figure [19] shows the evolution of a MBHB with q = 
1/81, Mi = 1, and crioo = 1 according to our hybrid 
model. The unit of time, Poj is the initial binary or- 
bital period at a = do, equal to (400,6700) yr when 
7 = (2, 1.5). When stellar disruptions arc not taken into 
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Fig. 20. — Stellar disruption rates as a function of time for an 
evolving MBHB with q = 1/81, M7 = 1, <rioo = 1, and different 
combinations of the initial e,7 values of the system. Each dotted 
curve is derived from 10 4 test scattering experiments in which the 
MBHB is evolved according to the rates give n by the hybrid model 
(see the s olid lines in the right panels of Fig. I19tl . Other line styles 
as in Fig. 1191 



account, the orbital semimajor axis shrinks by a factor 
of 10 on a timescale of 500Po before the binary stalls: at 
the same time, the eccentricity increases significantly to 
e ~ 0.5 — 1, depending on the initial value and of the 
parameter 7. When compared w ith the central panel of 
Figure 7 in ISesana et al.l (|2008|) , the results of the two 
integration schemes appear to be in excellent agreement. 
The inclusion of stellar disruptions causes the binary to 
stall at slightly larger a and higher e. The increase in 
the stalling radius is caused by the partial suppression 
of energetic ejections in favor of tidal disruptions. The 
larger eccent r icity increase can be explained as follow. 
ISesana et al.l (120081 ) demonstrated that stars with a* < a 
tend to drive the binary toward circularization, while 
stars with a* > a tend to increase its eccentricity. Since 
the former are the most susceptible to tidal disruptions, 
and disrupted stars do not exchange energy and angular 
momentum with the binary on the average, the relative 
contribution of stars with a* > a is larger in the pres- 
ence of tidal disruptions, pushing the binary eccentricity 
to higher values. In a realistic situation, the shrinking of 
the binary would continue at t > 500Po because of loss- 
cone diffusion processes and gravitational wave emission, 
which are not considered in this study. 

Figure HD] shows the stellar-disruption rates predicted 
by the hybrid model (solid and dashed lines) together 
with those derived by the test scattering experiments 
with the MBHB evolved by hand (dotted lines). Dur- 
ing the first 500Po, the rate remains constant, at a level 
that is comparable to the peak value calculated for a 
stationary binary. The duration of the plateau, how- 
ever, is longer in the case of a decaying pair, as new 
stars are continuously added to the time-varying inter- 
action loss cone. At t > 500Po> the evolution time of the 
MBHB exceeds the tidal disruption timescale, strongly- 
interacting stars get depleted, and the consumption rate 
drops sharply. The figure also shows that the peak dis- 
ruption rate is not sensitive to the binary eccentricity but 
depends on 7 according to the scaling relation in equa- 
tion (|3"7) . After 10 8 yr, the total number of disrupted 
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stars is (6.5 x 10 4 ,2.3 x 10 4 ) for 7 = (2,1.5). The dis- 
ruption rate during the plateau phase remains constant 
even though equation (|37[) predicts an increase oc ag/a 
(assuming an isothermal cusp). Equation (|37|) was de- 
rived, however, assuming a stationary binary at different 
orbital separations in an unperturbed stellar profile. In 
our hybrid model, the binary shrinks while depleting the 
stellar cusp, and many of the stars available for disrup- 
tion at, say, a = 0.1 ao are actually ejected or disrupted 
during the evolution of the pair from from ao to O.lao, 
leveling off the disruption rate. Given that: (1) the dis- 
ruption rate in the plateau phase remains constant and 
is consistent with the predictions of the Kozai mecha- 
nism even for evolving binaries; and (2) the duration of 
the plateau is of the order of the binary decay timescale, 
t d oc g~ 3 / 2 P (|Sesana et al.ll2008t ). the total number of 
disrupted stars can be scaled according to 

N. oc t d N* oc (3 - 7 )- 1 / 2 a( 2 -^/( 6 - 2 ^M 1 2/3 a, 

OC (3 - 7 )-V2 g (2- 7 )/(6-2 7 ) M ll/12 ) (47) 

where we used equation (|TT|) for Po and the Mi — <r* 
relation in the second proportionality. According to the 
above equation, for q = 1/81, should drop by a fac- 
tor of 2.5 as 7 varies from 2 to 1.5, consistent with the 
numbers derived from our hybrid model. Also, as long 
as 7 > 1.5, iV* should not be very sensitive to the binary 
mass ratio q. We stress that these scalings are derived 
from the no-GR/no-cusp experiments, and their validity 
is limited to binaries with q > 0.01. 

8. SUMMARY AND CONCLUSIONS 

In this paper, we have studied the tidal disruption rate 
in a system composed by a MBHB and a bound stellar 
cusp. We have carried out numerical scattering experi- 
ments for a detailed investigation of the mechanisms re- 
sponsible for the repopulation of the tidal loss cone, and 
developed a hybrid model to self-consistcntly solve for 
the evolutions of the binary, the depletion of the stellar 
cusp, and the stellar consumption rate. Our main results 
can be summarized as follows: 

1. For unequal binaries (q < 0.1), the tidal disrup- 
tion cross section for bound stars, which quanti- 
fies the probability of stellar disruption, is three 
orders of magnitude larger than the cross section 
for a single MBH fed by two-body relation. Two 
mechanisms contribute to such enhancement, the 
Kozai secular effect and chaotic resonant interac- 
tions. When the eccentricity of the MBHB is small, 
stars inside the Kozai wedge repopulate the tidal 
loss cone on the Kozai timescale, while stars out- 
side the Kozai wedge but inside the interaction loss 
cone are scattered into the tidal loss cone at ran- 
dom times due to close interactions with the sec- 
ondary hole. When the eccentricity is large, chaotic 
loss-cone repopulation becomes dominant over the 
entire range of stellar semimajor axis a* > (1 — e)a. 

2. GR and cusp- induced precession quench the Kozai 
secular evolution of interacting stars, causing a sig- 
nificant suppression (by a factor of ~ 10) of the 
disruption rate for q < 0.01. Therefore, the opti- 
mal enhancement of the tidal disruption rate by a 



MBHB occurs for mass ratios 0.01 < q < 0.1. Note 
that even if suppressed by a factor of ~ 10, the 
tidal disruption rate for binaries with q < 0.01 is 
still two order of magnitude larger than that given 
by standard relaxation processes around a single 
MBH. 

3. If a MBHB with mass ratio q <C 1 does not evolve 
significantly during 1/q revolutions, tidal disrup- 
tions of bound stars could initially persist at a 
constant rate ("plateau phase") that is four dex 
higher than the typical rates predicted for single 
MBHs. After one Kozai timescale (evaluated at 
a* = a), the tidal loss cone is repopulated mainly 
by chaotic interaction, and the stellar disruption 
rate decreases with time. The majority of stars are 
disrupted during a post-plateau later phase. 

4. If a MBHB evolves significantly on a timescale of 
1 1 q revolution, the plateau phase of stellar disrup- 
tions may last longer than a Kozai timescale. Tidal 
disruptions of bound stars slow down the shrinking 
of the binary and speed up the growth of binary 
eccentricity. 

Our results indicate that, after the formation of an 
unequal-mass MBHB at the center of a dense stellar cusp, 
the tidal disruption rate may go through three distinct 
evolutionary phases. The first phase begins shortly af- 
ter the MBHs become bound, and is characterized by a 
disruption rate as high as 0.1 — 1 stars per year, result- 
ing from the three-bo dy interactions be tween the binary 
and the bound stars (jChen et al.ll2009| ). When the de- 
cay timescale of the MBHB becomes longer than the tidal 
disruption timcscales of stars with a» ~ a, a second phase 
starts, where cusp depletion from slingshot ejections and 
tidal disruptions c auses a sharp drop in the disruption 
rate. iChen et al.l ()2008l ) showed that, unless stellar re- 
laxation is far more efficient than two-body "collisions" , 
the tidal disruption rate in this phase is orders of mag- 
nitudes lower than typical for a single MBHs. A third 
phase begins if the MBHB shrinks to the gravitational 
wave regime and eventually coalesces. The tidal disrup- 
tion rate then gradually increases to the value typical for 
single MBHs, 10 ~ 5 — 10~ 4 yr -1 , within o ne stellar relax- 
ation timescale (|Merritt fc Wang~l 120051 ). The number 
of stars disrupted during phase I is about 10 4 — 10 5 for 
AI7 = 1 and q = 1/81. The number of stars disrupted in 
phases II and III depends on the efficiency of stellar re- 
laxation, but would not significantly exceed ~ 10 5 — 10 6 . 
If a galaxy formed, on the average, one unequal-mass 
MBHB following a minor merger in its lifetime, then 
the above numbers imply that in a sample of tidal flares 
from MBHs of ~ 10 7 M Q , about 10% of events would be 
associated to binaries. If a galaxy forms unequal-mass 
MBHBs multiple times during its lifetime, then the de- 
tection rate of tidal events from binaries in transient sur- 
veys may be higher. Given the very high rates, there is 
also the possibility to identify MBHBs in galaxies host- 
ing multiple tidal flares within a years-to-decades time 
span. Over the next decade, synoptic surveys are ex- 
pected to detect hundreds of tidal disruption candidates. 
A tidal flare associated to a MBHB is l ikely interrupte d 
within one orbital period of the binary ()Liu et al.ll2009ft , 
therefore is distinguishable from the flares produced by 
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single MBHs, as long as the orbital period of the binary 
is shorter than the duration of a transient survey. If 
MBHB-drivcn disruptions account for 10% of the total 
rate, then the prospects of identifying MBHBs through 
tidal flares arc promising. Because the predicted rates in 
the three phases are significantly different from one an- 
other, the average stellar disruption rate over the lifetime 
of a galaxy is sensitive to the infalling rate of secondary 
MBHs and the relative duration of each phase. A com- 
parison tK3twej3nJdie<^ of tidal 
events (jDonlev et al.ll2002l : iGezari et aTll2008D and those 
predicted during the three phases may then shed light on 
the abundance and dynamical evolution of MBHBs. 
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